Theory of ultrafast nonequilibrium dynamics in d-wave superconductors 
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We use density-matrix theory to calculate the ultrafast dynamics of unconventional superconduc- 
tors from a microscopic viewpoint. We calculate the time evolution of the optical conductivity as 
well as pump-probe spectra for a d-wave order parameter. Three regimes can be distinguished in 
the spectra. The Drude response at low photon energies is the only one of those which has been 
measured experimentally so far. At higher energies, we predict two more regimes: the pair-breaking 
peak, which is reduced as Cooper-pairs are broken up by the exciting pulse; and a suppression above 
the pair-breaking peak due to nonequilibrium quasiparticles. Furthermore, we consider the influence 
of the electron-phonon coupling, and derive rate equations which have been widely used so far. 



PACS numbers: 74.25.Gz,74.72.-h,74.20.Rp,74.25.kc 

Introduction - In recent years, numerous studies of the 
nonequilibrium dynamics of carriers in superconductors 
have been performed using femtosecond time-resolved 
spectroscopy 0, 0, 0, 0, 0, H 0] ■ In a typical experiment, 
the sample is excited with an intense fs laser pulse (pump 
pulse), and after a delay time At, spectra are measured 
using a second, less intense, laser pulse (probe pulse). As 
the nature of the interactions between quasiparticles in 
the high-T c cuprates is still under debate [8J , it is interest- 
ing to directly observe the characteristic dynamics of con- 
densate depletion and Cooper-pair recombination. This 
can be done with real-time optical techniques. In the 
high-T c superconductor Bi2Sr2CaCu20s+a (BSCCO), for 
example, relaxation times of about 50 ps have been mea- 
sured the observed decay is two-component (biexpo- 
nential). 

Theoretical attempts to model these experiments have, 
on the one hand, used quasi-equilibrium models (so- 
called /i*, T* models, [§]) to describe the state excited 
by the pump pulse. On the other hand, rate equation 
approaches based on the phenomenological Rothwarf- 
Taylor model [1(3] have been used [H, 0, 0] to describe the 
recovery dynamics of the superconducting state. It is as- 
sumed that the dynamics are governed by the creation of 
high-energy phonons due to Cooper-pair recombination 
and subsequent phonon decay. So far, there has been no 
attempt to describe the excitation and relaxation dynam- 
ics on equal footing. As well, no microscopic description 
of the related time dynamics is available. 

In this Letter, we present a theory which can describe 
the femtosecond excitation and relaxation processes from 
a microscopic viewpoint. In particular, we consider high- 
Tc cuprates, using a realistic band structure and consider- 
ing coupling to two important phonon modes (breathing 
and buckling modes, which are strongly coupled to the 
superconducting CuC-2 planes [HI]). We employ the ap- 
proach of density-matrix theory, which has been used to 



some extent to describe ultrafast dynamics in semicon- 
ductors (see e.g. Ref. [H, O, [3]). 

Theory - We start from a Hamiltonian H = H sc + 
-fffieid + -ffphon, where H sc describes the superconducting 
state, -f/fieid gives the interaction with the classical elec- 
tromagnetic field, and -ffphon models the bare phonons 
and their interaction with the electrons. Explicitly we 
write 

H sc = ^(e k -M)4 s c ks + ( A * c tt c \l + h - c ) > (!) 
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In Eq. |T]), £k is a tight-binding band structure as mea- 
sured by Kordyuk et al. [1 51 ] , [i is the chemical potential, 
and Ak = Ao(cosfc K — cos/c y )/2 denotes a d-wave order 
parameter with Ao = 30 meV. c + and c are the electronic 
creation and annihilation operators, respectively. The j 
index counts the different phonon modes. In Eq. A q 
denotes the Fourier component of the vector potential the 
superconductor interacts with. It includes both the pump 
and probe fields; as the interaction with the pump field 
is nonlinear, the quadratic terms in A are needed. H p h on 
in Eq. ([3|) includes the bare phonons, having the disper- 
sion u)qj, and the electron-phonon interaction, described 
by the coupling matrix elements g p kj S - b + and b are the 
creation and annihilation operators for the phonons. We 
consider the important breathing and buckling phonon 
modes (l6|. These two modes are thought to be most 
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At= 20 ps 




FIG. 1: Plot of the excitation process. Starting with an equi- 
librium quasiparticle distribution (a^Ok) before the pump 
pulse (left panel), a nonequilibrium distribution is excited 
(middle panel), which then relaxes back into equilibrium 
(right panel). The inset shows the exiting pulse, a 50 fs Gaus- 
sian. The k-vectors lie in the 2D Cu02 plane; the plot shows a 
cut along the antinodal direction, from (0, 7r) to (n,n), cross- 
ing the Fermi level at &f. The temperature is 4 K. 



strongly coupled to the superconducting state, and thus 
the most relevant for scattering processes which can lead 
to relaxation of exited quasiparticles. 

We first perform a Bogoliubov transformation a k = 
uuc^ — WkC-ki, /?k = u kC-kj.+Wk c £f . Within the Heisen- 
berg picture we calculate equations of motion for the Bo- 
goliubov quasiparticle densities (a£ ak 2 )(t), (/3 k /3k 2 )(i), 
which correspond to the excited states of a supercon- 
ductor, and the anomalous expectation values {o£ffi£), 
which correspond to the condensate of Cooper-pairs. The 
current density is then given by 
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The first and last terms include Bogoliubov quasipar- 
ticle densities, thus describing the contribution of the 
normal part in a two-fluid-model. The second term, in- 
cluding anomalous expectation values, describes the con- 
densate response. Both Bogoliubov quasiparticle den- 
sities and anomalous expectation values can be calcu- 
lated for a given delay time At. As the probe field 
-E-probe = — *Wj4 pro be is known, the optical conductivity a 
can be calculated via (j)(q,w) — — iuja(q, u))A pT0 be(<l, w). 
Only the q-independent conductivity cr(q — » 0, w) will be 
considered. 

Equations of motion - In order to calculate a, the 
equations of motion for the Bogoliubov quasiparticle dis- 
tributions and anomalous expectation values have to be 
solved. They both couple to phonon-assisted quantities 



e.g. (aj 1+ri awo(blL r v j + b a j)). We now use second-order 
cluster expansion O, [ijj and calculate equations of mo- 
tion for the phonon-assisted quantities, which couple to 
4-point quantities such as (a£ + a^ak+qais). At this 
point, the hierarchy is broken down by factorizing the 
4-point quantities. The phonons are assumed to remain 
equilibrated (bath approximation, (bZjbqj) — ► Uqj with 
the Bosc distribution nqj) while the quasiparticles are 
excited and relax. The equations for the phonon-assisted 
quantities can then be solved, giving rise to a system 
of integro-differential equations. For example, the equa- 
tion for the Bogoliubov quasiparticle occupation (a^evk) 



reads: 



1 P 

d t (a+a k ) = --k • A q Afkq (K/?k) - 



+ E 



qj 



ds 



(1 + n qj )e ,;( " k +i^ k+ ^ )s 



x i kq UkMk+q(ak Q; k>(i - s)(l - (a£ +<l ak+q)(t - s)) 

J(u> k+ q— Wk-Wq^S (5) 



X 



i kq UkMk+q(aJ + a k+q )(< - s)(l - (a£a k )(t - s)) 



x Mk q Mk+qWk(/3k+ q /3k+q)(i - s)(a£ak)(t - s) 



With L k q = Uk+qUk + Vk+q^k, M kq = Uk + q«k 

being the relevant matrix elements, and to p 
where E p = 



«k+q«k 

= Ep/h, 

p + A p is the Bogoliubov quasiparticle 

dispersion. There are 4 equations for the 4 expectation 
values appearing in Eq. ([4]), all with a similar struc- 
ture. The full system is published elsewhere [l8j]. On 
this level, the equations are similar to the ones obtained 
within the Keldysh formalism, with the difference that 
here the nonequilibrium distributions (a^a k )(t — s) with 
their full time-dependences contribute. 

By using the Markovian approximation (lij . the inte- 
grals can be solved and one finds, for example, 



A q M kq ((a k /? k ) 
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= (1 + nqj )w k +qUkivkq<5(Wk+q _ W k + Wqj)i 
= Hqj'Uk+qWkikq^C^k+q — ^k — Wqj), 

= Uk+ q WkMkq<5(wk+q - ^k + Wqj). This is 



with 

r (l) 

i kqj 

r (2) 
i kqj 
r (3) 
1 kqj 

Boltzmann-type equation describing both in- and out- 
scattering with phonons (r^L, rj^l.) and Cooper-pair 

recombination (r kt ^-) processes. Finally, numerical so- 
lution yields the Bogoliubov quasiparticle distributions 
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and anomalous expectation values, and thus the optical 
conductivity. 

Results - Exciting the initial Bogoliubov quasiparticle 
distribution (a+a k ) = f(E k ) = (1 + exp (^k/fesT))" 1 
with a fs pump pulse, a nonequilibrium distribution is 
created, as shown in Fig. [TJ The biggest changes are 
around the Fermi energy, and the distribution is clearly 
non-thermal - quasiparticle weight is rearranged and con- 
sequently the condensate is also in a non-thermal state. 
Because of scattering with phonons, this nonequilibrium 
distribution can subsequently relax back into an equili- 
brated one. 

The probe conductivity after the pump and pump- 
probe spectra are obtained using the calculated Bogoli- 
ubov quasiparticle distributions as shown in Fig. [2] We 
can identify three regimes: the low-energy part (I) shows 
the Drude response, i.e. the response of the normal 
part in a two-fluid-modcl. The low-frequency power laws 
for d-wave superconductors are still obeyed after excita- 
tion. As Cooper-pairs are broken up by the pump pulse, 
thus generating Bogoliubov quasiparticlcs, the Drude re- 
sponse gets stronger. At higher energies ~ 2Ao one finds 
the pair- breaking peak (II). It gets shifted after pump- 
ing, as the superconducting state is depleted and Cooper- 
Pairs are broken up. Above the pair-breaking peak (re- 
gion III), the absorption, a ~ ai/uj, is suppressed. In 
an absorption process, Cooper-pairs have to be broken 
up, and the generated quasiparticles have to have empty 
states above 2A to be excited into. As a large number of 
quasiparticles are already excited due to the pump pro- 
cess, there are less states available than at equilibrium, 
which decreases the absorption. So far, only the Drude 
response part (I) has been measured experimentally 
In principle, however, the regimes II and III could be 
measured in THz pump-THz probe experiments. 

The enhancement of the Drude response and the shift 
of the pair-breaking peak are also found in a T* model, 
where the excited quasiparticle distribution is assumed to 
be an equilibrium distribution with an effective temper- 
ature T* Q. However, the suppression above 2 An. is not 
found within a T* model (see Fig. [3]). It is a nonequilib- 
rium effect - simply speaking, enhancing the temperature 
does not create enough Bogoliubov quasiparticles to fill 
a large number of states above 2Ao. 

Apart from pump-probe spectra, one can also look 
at the time evolution of the optical conductivity, which 
yields additional information about the recovery dynam- 
ics of the superconducting state. Fig. |4]shows the change 
in the conductivity Act = a(At) — o~o, where cto is the 
equilibrium conductivity (without a pump pulse). Act 
initially rises rapidly, as nonequilibrium quasiparticles 
are created and the superconducting state is depleted. 
After pumping, it decays. The overall timescale of this 
decay is given by the electron-phonon coupling, and thus 
faster for the breathing mode which is more strongly cou- 
pled to the electronic states. The decay is biexponential, 
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FIG. 2: Conductivity spectra for buckling (upper panel) and 
breathing (lower panel) modes. The equilibrium spectrum 
(without pump pulse) is shown along with spectra for different 
delay times. The inset shows the change in the real (red) and 
imaginary (blue) part of the conductivity, for At — 0.5 ps. It 
is Act = a(At) — <7n with the equilibrium conductivity <rn. 
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FIG. 3: Comparison of a calculated spectrum (At = 0) with a 
spectrum calculated using the effective T* model. T* is cho- 
sen in order to fit the calculated position of the pair-breaking 
peak. The inset shows the dependence of the T* model gap 
A* on the pump intensity. 



with the two timescales corresponding to quasiparticle- 
phonon scattering and Cooper-pair recombination. 

Derivation of rate equation approaches - Our micro- 
scopic approach can be used to derive a system of rate 
equations, which has been introduced by Kabanov et al. 
[l| to describe the combined dynamics of the excited 
quasiparticles and high-frequency phonons. The system 
is given by 



N 



I +r)N- Rn 2 

N n 2 
--Jv-V^+R^-l{N-N T ). 



(7) 



n, N are the numbers of exited quasiparticles and 
phonons, respectively, ?y, R are rates denoting pair- 
breaking and Cooper-pair recombination, and In, Jq are 
the initial changes in n and N. Nt is the equilibrium 
phonon number, and 7 describes phonon decay. 

So far, we have only considered phonons within the 
bath approximation, where they remain in equilibrium. 



4 




-80 ' L 1 1 1 1 -200 1 1 1 1 1 1 

20 40 60 80 100 2 4 6 8 10 

delay time A t [ps] delay time A t [ps] 



FIG. 4: Time evolution of the change in real (<n) and imag- 
inary (0-2) part of the optical conductivity for the buckling 
(left) and breathing (right) phonon modes. After an initial 
rise due to the pump process, a two-component decay of the 
conductivity follows. In region I scattering and in region II 
recombination dominate, respectively. 



phonon scattering and Cooper-pair recombination pro- 
cesses. The relaxation times calculated for the buckling 
modes are compatible with experimental results We 
have compared our results with spectra calculated within 
the T* model, finding good agreement in the low-energy 
limit, but our inclusion of nonequilibrium effects yields 
deviations at higher energies. Furthermore, we derive 
the widely used rate equation approaches from our mi- 
croscopic formalism. Our method thus provides insight 
into the condensate dynamics of d-wave superconductors 
and includes earlier theoretical attempts to describe it. 
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Our approach can be easily generalized to include the 
nonequilibrium phonon distributions within the Marko- 
vian approximation. The only modification in Eq. (JSJ) is 
in fact that the phonon distributions n q j are then time- 
dependent. A Boltzmann-like equation can also be de- 
rived for them. With n = ^ k (( a k a k) + {PfcP*)) = 
2 Ek( a it a k),we can derive an equation for n by summing 
Eq. ([6]) over all k. As only phonon absorption processes, 
i.e. pair-breaking by phonons, are relevant for Eq. (7]), 
only the first (initial values) and the last two terms in 
© need to be considered. The first term gives an initial 
rate I = £ k [k • A q M kq «a k /3 k ) - (a£/J+»] . As- 
suming constant recombination and phonon absorption 
rates, rg. - and rg. = n^-f^ - n q f t 2 ), one 
directly gets the form of Eq. (7J). A similar calculation 
with the phonon distribution equation yields the second 
rate equation. Thus, our microscopic approach includes 
the rate equations approach in the limit of constant scat- 
tering rates rW. We can then write the rates R and rj 
as: 

R = r( 3 >, V = f< 2 ) ]T<a+a k ) (1 - (a+Ok» , (8) 
k 

where the rates I?W are k, q, j-averages of the original 
rates, i.e. rW=E kqj -r k ^.. 

Conclusions - We have utilized density-matrix the- 
ory to calculate the ultrafast dynamics of high-T c su- 
perconductors. Our novel microscopic description of 
the optical excitation includes both the depletion of 
the superconducting condensate, as well as relaxation 
of the excited quasiparticles and Cooper-pair recombi- 
nation due to electron-phonon scattering. Pump-probe 
spectra, showing nonequilibrium effects above 2A , have 
been calculated as well as the real-time dynamics, where 
we find a biexponential decay produced by quasiparticle- 
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